#include "Helper.h"
#include "Spline.h"

using std::cout;
using std::endl;

int main(int argc, char* argv[]){
  cout << "------------------Problem B--------------------" << endl;
  int N[] = {5,10,20,40,80};
  cout << "For complete condition:" << endl;
  std::vector<Spline<1,4,ppForm> > Bv1;
  for (int i = 0 ; i < 5 ; i++){
     auto pvB = DataGenerator(Bfunc,-1,1,N[i]);
     Spline<1,4,ppForm> Sp = fitCurve<1,4>(pvB,complete,25.0/338,-25.0/338);
     Bv1.push_back(Sp);
  }
  Bplot(Bv1,"Output/B_complete.m");
  cout << "Run Output/B_complete.m by matlab to get the plot." <<
  endl;
  std::vector<double> err1(5);
  for (int i = 0 ; i < 5 ; i++)
    err1[i] = Berror(Bv1[i],N[i]);
  for (int i = 0 ; i < 4 ; i++){
    cout << "N =  " << N[i]+1 << " :\t error = " << err1[i] << ",\
  \t convergence rate = " << log(err1[i]/err1[i+1])/log(2) << endl;
  }
  cout << "N =  " << N[4]+1 << " : \t error = " << err1[4] << endl;
  cout << "For not-a-knot condition:" << endl;
  std::vector<Spline<1,4,ppForm> > Bv2;
  for (int i = 0 ; i < 5 ; i++){
     auto pvB = DataGenerator(Bfunc,-1,1,N[i]);
     Spline<1,4,ppForm> Sp = fitCurve<1,4>(pvB,notAknot);
     Bv2.push_back(Sp);
  }
  Bplot(Bv2,"Output/B_notAknot.m");
  cout << "Run Output/B_notAknot.m by matlab to get the plot." <<
  endl;
  std::vector<double> err2(5);
  for (int i = 0 ; i < 5 ; i++)
    err2[i] = Berror(Bv2[i],N[i]);
  for (int i = 0 ; i < 4 ; i++){
    cout << "N =  " << N[i]+1 << " :\t error = " << err2[i] << ",\
  \t convergence rate = " << log(err2[i]/err2[i+1])/log(2) << endl;
  }
  cout << "N =  " << N[4]+1 << " : \t error = " << err2[4] << endl; 
  cout << "------------------Problem C--------------------" << endl;
  cout << "For Corollary 4.58 case:" << endl;
  auto Pc1 = Cgenerator1();
  Spline<1,4,cardinalB> cSp1 = fitCurve<4>(-8,8,Pc1.first,Pc1.second);
  Cplot(cSp1,"Output/C_Bspline1.m");
  cout << "Run Output/C_Bspline1.m by matlab to get the plot." <<
  endl;
  cout << "For Corollary 4.59 case:" << endl;
  auto Pc2 = Cgenerator2();
  Spline<1,3,cardinalB> cSp2 = fitCurve<3>(-7,7,Pc2.first,Pc2.second);
  Cplot(cSp2,"Output/C_Bspline2.m");
  cout << "Run Output/C_Bspline2.m by matlab to get the plot." <<
  endl;
  cout << "------------------Problem D--------------------" << endl;
  double site[] = {-3.5,-3,-0.5,0,0.5,3,3.5};
  cout << "\t \t C_Bspline1 \t C_Bspline2" << endl;
  for (int i = 0 ; i < 7 ; i++){
    double tmperr1 = Derror<4>(cSp1,site[i]);
    double tmperr2 = Derror<3>(cSp2,site[i]);
    cout <<"Error at " << site[i] << " :" << "\t " << tmperr1;
    if (i == 3)
      cout << "\t";
    cout << "\t " << tmperr2 << endl;
  }
  cout << "We find that the error of C_Bspline1 reaches machine ";
  cout << "precision at -3, 0 and 3, while C_Bspline2 does the same at ";  
  cout << "-3.5, -0.5, 0.5 and 3.5. That's because these points are ";
  cout << "interpolation points for two splines respectively." <<
  endl;
  cout << "Notice that C_Bspline2 performes badly near 0, "; 
  cout << "therefore C_Bspline1 is more accurate." << endl;
  cout << "------------------Problem E--------------------" << endl;
  int num[] = {10,40,160};
  const Vec<double,2> Dbegin({0,1});
  const Vec<double,2> Dend({0,-1});
  std::vector<std::pair<Spline<2,4,ppForm>,Spline<2,4,ppForm> > > Evpsp;
  for (int i = 0 ; i < 3 ; i++){
    auto info = Egenerator(num[i]);
    Spline<2,4,ppForm> Esp1 =
  fitCurve<2,4>(info.first,complete,Dbegin,Dend);
    Spline<2,4,ppForm> Esp2 =
  fitCurve<2,4>(info.second,complete,Dbegin,Dend);
    Evpsp.push_back(std::make_pair(Esp1,Esp2));
  }
  Eplot(Evpsp,"Output/E_heart.m",400);
  cout << "Run Output/E_heart.m by matlab to get the plot." <<
  endl;
  cout << "------------------Problem F--------------------" << endl;
  InterpConditions<double,2> IC1("Input/Curve1");
  InterpConditions<double,2> IC2("Input/Curve2");
  InterpConditions<double,2> IC3("Input/Curve3");
  std::vector<Spline<1,4,ppForm> > Fvsp;
  Spline<1,4,ppForm> Fsp1 = interpolate<4,2,ppForm>(IC1,complete);
  Spline<1,4,ppForm> Fsp2 = interpolate<4,2,ppForm>(IC2,complete);
  Spline<1,4,ppForm> Fsp3 = interpolate<4,2,ppForm>(IC3,complete);
  Fvsp.push_back(Fsp1);
  Fvsp.push_back(Fsp2);
  Fvsp.push_back(Fsp3);
  cout << "(In the following display, all polynomials take the left ";
  cout << "end of the interval as the point that t = 0, while t here ";
  cout << "denote x)" << endl;
  cout << "Curve 1 :" << endl;
  Fsp1.show();
  cout << "Curve 2 :" << endl;
  Fsp2.show();
  cout << "Curve 3 :" << endl;
  Fsp3.show();
  Fplot(Fvsp,"Output/F_spline.m");
  cout << "Run Output/F_spline.m by matlab to get the plot." <<
  endl;
  cout << "------------------Problem G--------------------" << endl;
  Gplot<1>("Output/G_table1.m",80);
  Gplot<2>("Output/G_table2.m",160);
  cout << "For n = 1, Run Output/G_table1.m by matlab to get ";
  cout << "the picture." << endl;
  cout << "For n = 2, Run Output/G_table2.m by matlab to get ";
  cout << "the picture." << endl;
}
